Descriptions of dataset¶
- tas_df: temperature observations of 9 locations
- pr_df: precipitation observations of 9 locations
- flow_df: the response variable that we need to fit and predict
- Time: from 1958-01-04 to 2015-12-31 daily data of 21778 days
In [41]:
df.info()
In [42]:
df.describe()
Out[42]:
In [79]:
# timeseries normalized?
# strong annual cycle embedded >> need to remove before any fitting
# any trends?
tas_df.plot(subplots=True, layout=(9, 1), figsize=(8, 25), sharex=True) # certainly not normalized
sns.despine()#)
plt.show()
In [88]:
tas_df[tas_df['tas_0']<-30]
# The cold days seem quite universal across locations, for example, 1987-01-10 to 1987-01-11, and 1978-12-31 to 1979-01-01
# The cold days seem quite universal across locations, for example, 1987-01-10 to 1987-01-11, and 1978-12-31 to 1979-01-01
Out[88]:
In [103]:
#label=['Loc1','Loc2','Loc3','Loc4','Loc5','Loc6','Loc7','Loc8','Loc9']
ax=tas_df.plot(kind='box')
#ax.set_xticklabels(label,ha='center')
#ax = sns.boxplot(data=tas_df)
plt.title('Air_temperature spreads')
plt.ylabel('degC')
plt.show()
In [94]:
tas_df.corr() ## very high correlation because they share a similar annual cycle,
## The air_temperatures between locations appear to correlate pretty well.
Out[94]:
In [99]:
sns.pairplot(tas_df)
Out[99]:
In [13]:
# Remove seasonal cycle then estimate correlation
climatology=tas_df.groupby(tas_df.index.month).mean()
climatology.plot() ## seasonal cycle
Out[13]:
In [14]:
anomaly = lambda x: (x - x.mean())
deseasoned_tas=tas_df.groupby(tas_df.index.month).transform(anomaly)
In [27]:
tas_df['tas_0']['1960':'1965'].plot()
Out[27]:
In [26]:
deseasoned_tas['tas_0']['1960':'1965'].plot()
Out[26]:
In [28]:
sns.pairplot(deseasoned_tas)
Out[28]:
In [29]:
deseasoned_tas.corr()
Out[29]:
Linear trends in time¶
In [538]:
deseasoned_df.apply(linTrend)*365*10
Out[538]:
Precipitation¶
In [78]:
pr_df.plot(subplots=True, layout=(9, 1), figsize=(8, 25), sharex=True) # certainly not normalized
sns.despine()#)
plt.show()
In [127]:
# Remove seasonal cycle then estimate correlation
climatology_pr=pr_df.groupby(pr_df.index.month).mean()
climatology_pr.plot() ## seasonal cycle
# Most precipitation occurs in Summer, peaks in July, August, but remain high till November
# Lowest precipitation in February till April
Out[127]:
In [129]:
anomaly2 = lambda x: (x - x.mean())
deseasoned_pr=pr_df.groupby(pr_df.index.month).transform(anomaly2)
sns.pairplot(deseasoned_pr)
Out[129]:
In [132]:
deseasoned_pr.apply(linTrend)*365*10 #also a positive trends in precipitation in all locations
Out[132]:
Precipitation Summary¶
- Location 6-9 have high correlation, especially the high precip events in the first decade
- All locations have a high precip peak on 1994-08-14
Correlation of deseasoned data¶
In [141]:
sns.pairplot(deseasoned_df)
Out[141]:
In [157]:
for TS1, TS2 in zip(pr_df.columns, tas_df.columns):
print(deseasoned_pr[TS1].corr(deseasoned_tas[TS2]))
# After removing seasonal cycles, the correlations between rainfall and temperature at each location is very low.
River Flow¶
- Trend of the flow
- Seasonal cycle of the flow (should be one close to Precip)
- Find the location with greatest correlation with the Flow
In [20]:
# Looks very much like my Agulhas leakage timeseries.
sns.set(font_scale=2)
sns.set_style(style='white')
fig=plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
# the size of A4 paper
flow_df.plot(ax=ax) #
plt.ylim([0,90])
sns.despine()#)
plt.show()
In [360]:
# Remove seasonal cycle then estimate correlation
climatology_fl=flow_df.groupby(flow_df.index.month).mean()
climatology_fl.plot() ## seasonal cycle
# Highest river flow events in May.
Out[360]:
Spectral analysis of flow¶
In [560]:
from mtspec.multitaper import mtspec
PSD,f=mtspec(flow_df, 1, 4, number_of_tapers=7)
In [562]:
sns.set_style('whitegrid')
plt.semilogy(f,PSD)
plt.xlim([0,0.02]) # 0.5
plt.xlabel('Frequency(cycles/day)')
Out[562]:
Summary of the findings:¶
- Strong seasonal signals in precipitation, temperature in all locations, as well as the flow
- But the Tas and Pr peak in Summer (June to August), yet the Flow peaks in May
- All have positive linear trends in times, but not sure if such trends are significant (0.3deg/decade, 0.04mm/decade, 0.15m^3/decade).
- Low correlation between precipitation and temperature at same locations
- But high correlation among Temperature/precipittions among different locations, respectively.
- Correlation between tas and flow ~0.3, precip and flow ~0.2
Machine learning model > use 18 features to predict the river flow¶
- Try simple regression first (multivariate)
- Remove co-linearity between precip/temp
- Normalize/remove seasonal cycle before fitting
- Feature engineering, regularization
- Evaluate model performance
Naive approach, full dataframe¶
In [370]:
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(noflow, df['flow'], test_size=0.3, random_state=0)
from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train, y_train)
y_pred = regressor.predict(X_test)
print('Linear Regression R squared": %.4f' % regressor.score(X_test, y_test))
How to improve?¶
- Scaling the features
- cross-validation
- Regularization
- hyperparameter tuning use GridSearchCV
- remove seasonal cycle from the dataset
- Try different algorithms?
In [443]:
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
X_train, X_test, y_train, y_test = train_test_split(noflow, df['flow'], test_size=0.3, random_state=0)
steps = [('scaler', StandardScaler()),
('reg', LinearRegression())]
# Create the pipeline: pipeline
pipeline = Pipeline(steps)
pipeline.fit(X_train, y_train)
# Show scores:
from sklearn.model_selection import cross_val_score
cv_results=cross_val_score(pipeline,X_train, y_train,cv=5)
print(cv_results)
In [433]:
# Regularization Ridge Regression
from sklearn.linear_model import Ridge
steps2 = [('scaler', StandardScaler()),
('Ridge', Ridge(alpha=0.1,normalize=True))]
# Create the pipeline: pipeline
pipeline2 = Pipeline(steps2)
pipeline2.fit(X_train, y_train)
# Show scores:
cv_results=cross_val_score(pipeline2,X_train, y_train,cv=5)
print(cv_results)
In [444]:
# Regularization Lasso Regression
from sklearn.linear_model import Lasso
steps3 = [('scaler', StandardScaler()),
('Lasso', Lasso(alpha=0.1))]
# Create the pipeline: pipeline
pipeline3 = Pipeline(steps3)
pipeline3.fit(X_train, y_train)
# Show scores:
cv_results=cross_val_score(pipeline3,X_train, y_train,cv=5)
print(cv_results)
# Most important features:
Lasso_coef=pipeline3.named_steps['Lasso'].coef_
names=noflow.columns
_=plt.plot(range(len(names)),Lasso_coef)
_=plt.xticks(range(len(names)),names,rotation=60)
_=plt.ylabel('coefficients')
plt.show()
# Seems tas8, pr3, pr7, pr8 most representative for flow prediction.
In [435]:
# Remove Seasonal cycle
anomaly = lambda x: (x - x.mean())
deseasoned_df=df.groupby(df.index.month).transform(anomaly)
deseasoned_noflow=deseasoned_df.drop(['flow'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(deseasoned_noflow, deseasoned_df['flow'], test_size=0.3, random_state=0)
pipeline3.fit(X_train, y_train)
# Show scores:
from sklearn.model_selection import cross_val_score
cv_results=cross_val_score(pipeline3,X_train, y_train,cv=5)
print(cv_results)
# For all three Linear Regression, the R^2 score reduces significantly after removing seasonal cycle
In [463]:
# Use GridSearchCV to tune hyperparameters
from sklearn.model_selection import GridSearchCV
param_grid={'Lasso__alpha':np.linspace(0.1,2,20)}
lasso_cv=GridSearchCV(pipeline3,param_grid,cv=5)
lasso_cv.fit(X_train,y_train)
lasso_cv.best_score_
Out[463]:
Random Forest¶
In [468]:
from sklearn.ensemble import RandomForestRegressor
X_train, X_test, y_train, y_test = train_test_split(noflow, df['flow'], test_size=0.3, random_state=0)
#X_train, X_test, y_train, y_test = train_test_split(deseasoned_noflow, deseasoned_df['flow'], test_size=0.3, random_state=0)
steps = [('scaler', StandardScaler()),
('forest', RandomForestRegressor(random_state=42))]
# Create the pipeline: pipeline
pipeline = Pipeline(steps)
pipeline.fit(X_train, y_train)
print('Random Forest R squared": %.4f' % pipeline.score(X_test, y_test))
Gradient Boosting¶
In [471]:
from sklearn import ensemble
from sklearn.ensemble import GradientBoostingRegressor
X_train, X_test, y_train, y_test = train_test_split(noflow, df['flow'], test_size=0.3, random_state=0)
#X_train, X_test, y_train, y_test = train_test_split(deseasoned_noflow, deseasoned_df['flow'], test_size=0.3, random_state=0)
steps = [('scaler', StandardScaler()),
('gb', GradientBoostingRegressor())]
# Create the pipeline: pipeline
pipeline = Pipeline(steps)
pipeline.fit(X_train, y_train)
print('Gradient Boost R squared": %.4f' % pipeline.score(X_test, y_test))
Different Approach: Classification of high-flow events rather than regression¶
- Labeling df['flow'] by prescribing a threshold
- Try KNN and Logistic
In [472]:
df['flow'].plot()
Out[472]:
In [475]:
df['flow'].describe(percentiles=[.25, .5, .75,.90,.95])
Out[475]:
In [507]:
# Labeling for flow > 20
df['flood'] = np.where(df['flow']>=20, 1,0)
df_labeld=df.drop(['flow','flood'], axis=1)
KNN¶
In [565]:
from sklearn.neighbors import KNeighborsClassifier
X_train, X_test, y_train, y_test = train_test_split(df_labeld, df['flood'], test_size=0.3, random_state=21)
knn=KNeighborsClassifier(n_neighbors=6)
knn.fit(X_train,y_train)
y_pred=knn.predict(X_test)
y_prob=knn.predict_proba(X_test)
print(knn.score(X_test,y_test))
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))
Logistic Classification¶
In [514]:
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix, classification_report
X_train, X_test, y_train, y_test = train_test_split(df_labeld, df['flood'], test_size=0.3, random_state=21)
# Create the classifier: logreg
logreg = LogisticRegression()
# Fit the classifier to the training data
logreg.fit(X_train,y_train)
# Predict the labels of the test set: y_pred
y_pred = logreg.predict(X_test)
# Compute and print the confusion matrix and classification report
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))
In [535]:
from sklearn.metrics import roc_curve
# Compute predicted probabilities: y_pred_prob
y_pred_prob = logreg.predict_proba(X_test)[:,1]
# Generate ROC curve values: fpr, tpr, thresholds
fpr, tpr, thresholds = roc_curve(y_test, y_pred_prob)
# Plot ROC curve
plt.plot([0, 1], [0, 1], 'k--')
plt.plot(fpr, tpr)
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
plt.show()
In [537]:
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV
# Create the hyperparameter grid
c_space = np.logspace(-5, 8, 15)
param_grid = {'C': c_space, 'penalty': ['l1', 'l2']}
# Instantiate the logistic regression classifier: logreg
logreg = LogisticRegression()
# Create train and test sets
X_train, X_test, y_train, y_test = train_test_split(df_labeld, df['flood'], test_size=0.3, random_state=21)
# Instantiate the GridSearchCV object: logreg_cv
logreg_cv = GridSearchCV(logreg,param_grid,cv=5)
# Fit it to the training data
logreg_cv.fit(X_train,y_train)
# Print the optimal parameters and best score
print("Tuned Logistic Regression Parameter: {}".format(logreg_cv.best_params_))
print("Tuned Logistic Regression Accuracy: {}".format(logreg_cv.best_score_))
Next steps¶
- Hyperparameter tuning
- Try other algorithms, e.g. Neural Networks
- Better way to split the data (much more non-flood days vs flood days)
- Not sure if the annual cycle should be considered.
- Feature selection, engineering > leave out some less representative locations